clear all
capture log close
global dir "directory"
log using "$dir/acamort/code/0Y_MJW_Roth.log", replace
/*******************************************************************************
PROGRAM:	09_MJW_Roth.do
PURPOSE:	Calculations from Roth (2019) for QJE R&R
DATE:		June 25, 2020
*******************************************************************************/
global indata "directory/indata/rdcoutput"
set seed 55739634

**read in estimates from RDC for main specification
import excel using "$indata/varcov.xlsx", clear first
mkmat Eventm5 Eventm4 Eventm3 Eventm2 Eventm1 Event1 Event2 Event3 Event4, mat(v) rown(A)
mat list v
mat v2 = vecdiag(v)'

import excel using "$indata/beta.xlsx", clear
mkmat B, mat(b) rown(A) nomissing
mat list b

**calculate slope with rejection probability of 80%
**as described on pg. 23, assume pre-trend coefficients drawn from 
**normal distribution where mean is a vector of means with an upward slope gamma,
**and the variance is the component of the covariance matrix corresponding with the
**pre-trends coefficients
mat v_pre=v[1..5,1..5]
mat b_pre=b[1..5,1]
mat b_0=[-5,-4,-3,-2,-1]
mat b_0=b_0'
clear
*estimate trend at 80% as -0.0003235
foreach s in -.0003230 -.0003231 -.0003232 -.0003233 -.0003234 -.0003235 ///
-.0003236 -.0003237 -.0003238 -.0003239   {
	clear
	mat s=b_0*`s'
	drawnorm b1 b2 b3 b4 b5, means(s) cov(v_pre) n(100000)
	forval x=1/5{
		local lb=b[`x',1]-1.96*sqrt(v2[`x',1])
		local ub=b[`x',1]+1.96*sqrt(v2[`x',1])
		gen reject`x'=1*(~inrange(b`x',`lb',`ub'))
	}
	egen totreject=rowtotal(reject1-reject5)
	gen anyreject=totreject>0
	su anyreject
	drop b1 b2 b3 b4 b5 anyreject totreject reject1-reject5
}
*estimate trend at 50% as -0.0001774
foreach s in -.0001770 -.0001771 -.0001772 -.0001773 -.0001774 -.0001775 ///
-.0001776 -.0001777 -.0001778 -.0001779  {
	clear
	mat s=b_0*`s'
	drawnorm b1 b2 b3 b4 b5, means(s) cov(v_pre) n(100000)
	forval x=1/5{
		local lb=b[`x',1]-1.96*sqrt(v2[`x',1])
		local ub=b[`x',1]+1.96*sqrt(v2[`x',1])
		gen reject`x'=1*(~inrange(b`x',`lb',`ub'))
	}
	egen totreject=rowtotal(reject1-reject5)
	gen anyreject=totreject>0
	su anyreject
	drop b1 b2 b3 b4 b5 anyreject totreject reject1-reject5
}

**calculate unconditional bias and conditional bias for trend at 80% power
**Roth says to use Proposition 3.1 on page 16 to do this
*unconditional bias
mat b_1=[1,2,3,4]
mat uncondbias=b_1'*-.0003235
mat list uncondbias

*conditional bias
clear
mat s=b_0*-.0003235
drawnorm b1 b2 b3 b4 b5, means(s) cov(v_pre) n(100000)
forval x=1/5{
	local lb=b[`x',1]-1.96*sqrt(v2[`x',1])
	local ub=b[`x',1]+1.96*sqrt(v2[`x',1])
	gen reject`x'=1*(~inrange(b`x',`lb',`ub'))
}
egen totreject=rowtotal(reject1-reject5)
gen anyreject=totreject>0
su anyreject
keep if anyreject==0
collapse b1 b2 b3 b4 b5
mkmat b1 b2 b3 b4 b5, mat(mean)
mat mean=mean'
mat list mean
mat list s
mat distort=mean-s
mat list distort
mat v_post_pre=v[6..9,1..5]
mat list v_post_pre
mat inv_v_pre=inv(v_pre)
mat list v_pre
mat list inv_v_pre
mat stdcov=v_post_pre*inv_v_pre
mat condbias=(stdcov*distort)+uncondbias

**Show bias for each treatment effect
mat list b
mat list condbias

**calculate unconditional bias and conditional bias for trend at 50% power
*unconditional bias
mat uncondbias=b_1'*-.0001774
mat list uncondbias
*conditional bias

clear
mat s=b_0*-.0001774
drawnorm b1 b2 b3 b4 b5, means(s) cov(v_pre) n(100000)
forval x=1/5{
	local lb=b[`x',1]-1.96*sqrt(v2[`x',1])
	local ub=b[`x',1]+1.96*sqrt(v2[`x',1])
	gen reject`x'=1*(~inrange(b`x',`lb',`ub'))
}
egen totreject=rowtotal(reject1-reject5)
gen anyreject=totreject>0
su anyreject
keep if anyreject==0
collapse b1 b2 b3 b4 b5
mkmat b1 b2 b3 b4 b5, mat(mean)
mat mean=mean'
mat list mean
mat list s
mat distort=mean-s
mat list distort
mat v_post_pre=v[6..9,1..5]
mat list v_post_pre
mat inv_v_pre=inv(v_pre)
mat list v_pre
mat list inv_v_pre
mat stdcov=v_post_pre*inv_v_pre
mat condbias=(stdcov*distort)+uncondbias

**Show bias for each treatment effect
mat list b
mat list condbias
